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LONG-TERM  GOALS 

The  long-term  goal  of  this  research  is  to  construct  the  Navy’s  next-generation  global  numerical 
weather  prediction  (NWP)  model  using  new  numerical  methods  specifically  designed  for  distributed- 
memory  and  vector  computers.  To  distinguish  it  from  the  current  Navy  global  atmospheric  model,  the 
navy  Operational  Global  Atmospheric  Prediction  System  (NOGAPS),  we  shall  refer  to  this  new  model 
as  the  Navy’s  Spectral  Element  Atmospheric  Model  (NSEAM).  To  take  full  advantage  of  distributed- 
memory  computers,  the  spherical  global  domain  of  NSEAM  is  partitioned  into  local  sub-domains,  or 
elements,  which  can  then  be  solved  independently  on  multiple  processors.  The  numerical  methods 
used  on  these  sub-domains  must  be  not  only  local  in  nature  but  also  high-order  accurate  and  highly 
efficient.  Thus,  NSEAM  is  being  constructed  so  that  it  is  as  accurate  as  the  current  spectral  model 
(NOGAPS);  is  more  efficient,  thereby  allowing  for  finer  resolution  forecasts;  and  is  geometrically 
more  flexible,  thereby  allowing  for  the  use  of  adaptive  or  telescoping  grids.  This  will  allow  for  better 
coupling  with  mesoscale  models  and  eventually  perhaps  allowing  for  a  modeling  paradigm  whereby 
the  global  and  mesoscale  models  are  either  unified  or  at  least  share  the  same  numerical  algorithms 
which  facilitates  their  maintenance. 

OBJECTIVES 

The  objective  of  this  project  is  to  construct  new  high-order  local  methods  for  the  Navy’s  next- 
generation  global  NWP  model.  The  high-order  accuracy  of  these  methods  will  ensure  that  the  new 
model  yields  the  same  accuracy  as  the  current  spherical  harmonics  model  while  the  locality  of  these 
methods  will  ensure  that  the  efficiency  of  the  model  increases  as  the  number  of  processors  increases. 

APPROACH 

To  meet  our  objectives  we  explore: 

1 .  spectral  element  (SE)  and  discontinuous  Galerkin  (DG)  methods  on  quadrilateral  and  triangular 
grids, 

2.  semi-implicit  (SI)  and  semi-Lagrangian  (SL)  time-integrators,  and 

3.  adjustments  to  the  physical  parameterization  and  non-reflecting  boundary  conditions  of  the  current 
NOGAPS  to  incorporate  it  into  NSEAM. 
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The  power  of  SE  and  DG  methods  is  that  they  are  high-order  accurate,  like  spherical  harmonics 
methods,  yet  are  completely  local  in  nature  -  meaning  that  the  equations  are  solved  independently 
within  each  individual  element  and  processor.  Furthermore,  high-order  methods  have  minimal 
dispersion  error  which  is  an  important  property  for  capturing  fine-scale  atmospheric  phenomena  (e.g., 
tropical  cyclones,  Kelvin  and  Rossby  waves).  In  addition,  SI  and  SL  methods  offer  vast  improvements 
in  efficiency  due  to  the  longer  time  steps  that  they  permit. 

After  validation  of  the  SE,  DG,  SI,  and  SL  discretization  on  barotropic  test  cases,  the  full  3D  adiabatic 
(dry  physics)  hydrostatic  primitive  equation  model  was  closely  scrutinized.  The  SI  time-integrator  of 
NSEAM  was  ported  to  the  message-passing  interface  (MPI).  At  this  point  direct  comparisons  with  the 
current  version  of  NOGAPS  were  made.  Superior  performance  of  NSEAM  over  NOGAPS  then 
justified  the  further  effort  into  taking  the  steps  to  make  NSEAM  the  next  operational  NWP  model.  In 
order  to  accomplish  the  final  objective,  the  following  steps  must  be  accomplished: 

1 .  the  NOGAPS  physical  parameterization  package  must  be  modified  and  then  coupled  with  the 
NSEAM  dynamical  core, 

2.  non-reflecting  boundary  conditions  need  to  be  added  to  NSEAM  in  order  to  avoid  spurious  wave 
reflections  emanating  from  the  orography,  and 

3.  development  of  the  tangent  linear  and  adjoint  models  of  NSEAM  to  be  used  for  both  data 
assimilation  and  sensitivity  studies. 

WORK  COMPLETED 

NSEAM  Scalability.  The  superior  performance  of  the  semi-implicit  NSEAM  were  published  in  the  July 
B  2005  issue  of  the  Quarterly  Journal  of  the  Royal  Meteorological  Society  (Ref.  [4])  which  showed 
that  at  T239  L30,  NSEAM  dynamics  is  60%  faster  than  NOGAPS  dynamics,  and  that  at  T498  L60 
NSEAM  scales  linearly  -  a  feat  unequaled  by  spherical  harmonics  and  most  newly-proposed  next- 
generation  models  (see  Fig.  1). 

NSEAM  Numerics.  The  DG  method  on  triangular  grids  was  developed  (see  Ref.  [2])  which  shows  that 
this  method  is  capable  of  handling  sharp  discontinuities  (such  as  shocks)  which  then  makes  it  the 
leading  candidate  for  resolving  all  types  of  mesoscale  phenomena.  In  addition,  a  new  semi-implicit 
Lagrangian  time-integrator  has  been  developed  (see  Ref.  [1]),  which  we  will  now  begin  testing  on 
NSEAM;  this  new  time-integrator  should  allow  larger  time-steps  to  be  used  which  will  result  in  a 
further  increase  in  efficiency. 

Non -reflecting  Upper  Boundary  Conditions.  Non-reflecting  upper  boundary  conditions  must  be 
implemented  within  NSEAM  in  order  to  avoid  vertically  propagating  gravity  waves  emanating  from 
the  orography  to  reflect  downward  from  the  model  top  and  then  contaminating  the  solution.  From  the 
various  choices  available,  it  was  determined  that  using  the  conventional  approach  of  using  sponge 
layers  which  applies  Rayleigh  friction  to  the  surface  layers  and  Newtonian  cooling  to  the  top  layers 
would  result  in  the  fastest  working  model.  This  sponge  has  now  been  added  to  NSEAM  and  testing  of 
the  sponge  layer  has  begun.  To  test  the  effectiveness  of  the  sponge  ,we  have  added  analytic  mountains 
in  order  to  simulate  real  waves.  To  increase  the  effectiveness  of  the  sponge,  we  have  increased  the 
model  layers  to  48  but  expect  to  further  increase  it  to  60  in  2006. 
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Physical  Parameterization.  A  majority  of  the  time  in  this  project  this  year  has  been  spent  on  coupling 
the  NSEAM  dynamical  core  with  the  current  NOGAPS  physical  parameterization  routines  (called 
DIABAT).  Various  interfaces  have  had  to  be  constructed  in  order  to  allow  both  of  these  quite  different 
models  to  work  in  unison.  For  example,  NSEAM  relies  on  Cartesian  coordinates  to  allow  for  any  grid 
to  be  used  whereas  DIABAT  assumes  spherical  coordinates.  Thus  a  very  accurate  way  of  mapping 
from  Cartesian  to  spherical  coordinates  had  to  be  developed.  This  transformation  is  not  as 
straightforward  as  one  might  think  due  to  the  singular  nature  of  trigonometric  functions  near  the  poles. 
With  this  transfonnation  developed,  NSEAM  uses  Cartesian  coordinates  for  the  dynamics,  maps  to 
spherical  coordinates  for  the  physics  (DIABAT),  and  then  maps  back  to  Cartesian  for  the  remainder  of 
the  operations. 

Another  challenge  with  including  DIABAT  within  NSEAM  is  that  the  logic  of  these  two  routines  is 
vastly  different.  For  example,  NSEAM  is  written  in  a  general  way  using  indirect  addressing  which 
allows  for  any  grid  geometry  to  be  used.  However,  DIABAT  is  written  in  a  logically  Cartesian  way 
where  there  is  an  “I”  and  “J”  direction;  this  is  called  direct  addressing.  Direct  addressing  is  typical  of 
structured  grid  codes  whereas  NSEAM  is  written  for  fully  unstructured  grids  and  thereby  has  no 
direction.  This  incompatibility  between  the  two  codes  has  required  a  careful  modification  of  some 
portions  of  the  DIABAT  code  which  we  describe  in  the  results  section. 

RESULTS 

NSEAM  Scalability.  In  Ref.  [4]  we  show  that  NSEAM  scales  extremely  well  for  increasing  processor 
count,  especially  as  the  resolution  is  increased.  Fig.  1  shows  that  at  T249  L30  and  using  a  time-step  of 
300  seconds,  NSEAM  achieves  over  60  simulation  days  per  wallclock  hour  for  386  processors. 
However,  it  is  evident  from  Fig.  la  that  the  perfonnance  of  NSEAM  T249  L30  is  beginning  to  peak. 
Increasing  the  resolution  to  T498  L60  with  a  time-step  of  150  seconds  increases  the  problem  size  by  a 
factor  of  16  (a  factor  of  4  in  the  horizontal,  a  factor  of  2  in  the  vertical,  and  another  factor  of  2  in  time); 
however,  Fig.  lb  shows  that  NSEAM  T498  L60  achieves  over  7  simulation  days  per  wallclock  hour 
for  386  processors.  Thus  the  perfonnance  has  only  decreased  by  a  factor  of  8  even  though  the  problem 
size  has  increased  by  a  factor  of  16.  Furthermore,  NSEAM  has  scaled  linearly  (perfectly)  at  T498  L60. 


Figure  1:  Scalability  of  NSEAM  on  an  IBM  SP 4  with  resolution  a)  T249  L30  and  b)  T498  L60. 
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NSEAM  Numerics.  To  show  the  advantages  of  using  the  new  Lagrangian  semi-implicit  time-integrator 
that  we  developed  in  Ref.  [1]  in  Fig.  2,  we  show  error  norms  as  a  function  of  Courant  number  which  in 
turn,  is  a  function  of  time-step;  thus  a  large  Courant  number  signifies  a  large  time-step.  Thus  in  these 
figures,  the  Eulerian  semi-implicit  time-integrator  developed  previously  in  Ref.  [4]  can  only  use 
Courant  numbers  less  than  3,  which  is  the  first  point  on  the  left  in  Figs.  2a  and  2b.  The  new  time- 
integrator,  which  we  call  HELSI  (Hybrid  Eulerian-Lagrangian  Semi-Implicit),  combines  the  semi- 
Lagrangian  (SL)  and  operator-integration-factor  splitting  (OIFS)  methods  into  a  unified  formulation. 
From  Fig.  2  one  can  see  that  both  methods  yield  very  good  error  norms  for  Courant  numbers  up  to  15; 
in  other  words,  we  have  been  able  to  increase  the  time-step  by  a  factor  of  5  from  the  Eulerian  semi- 
implicit  method  without  incurring  large  errors.  While  this  magnitude  of  increased  time-step  may  not 
be  realistic  for  more  complicated  flow  problems  (with  terrain  and  highly  nonlinear  effects)  even  a 
modest  time-step  increase  of  a  factor  of  2  is  sufficient  to  make  the  new  approach  worthwhile.  Tests  on 
more  complicated  flows  will  be  reported  on  in  an  upcoming  journal  article  currently  being  written. 


Figure  2:  Error  as  a  function  of  Courant  number  for  the  new  time-integrators  using  a)  the  semi- 
Lagrangian  (SL)  and  b)  the  operator-integration-factor  splitting  (OIFS)  methods. 


Physical  Parameterization.  To  elegantly  splice  DIABAT  and  NSEAM  it  has  been  critical  to  modify 
some  portions  of  the  DIABAT  code,  in  addition  to  constructing  interface  routines.  The  difficulty 
stems  from  the  vastly  different  fonnulations  of  structured  and  fully  unstructured  grid  codes. 
Specifically,  because  DIABAT  assumes  a  latitude-longitude  grid,  which  is  a  logically  Cartesian  grid,  it 
expects  the  data  to  be  ordered  in  an  I  (longitude)  and  J  (latitude)  direction.  Thus  the  point  1-1  and  1+1 
of  any  variable  F  will  be  the  left  and  right  neighbors  of  the  point  I;  this  is  an  example  of  direct 
addressing.  In  contrast  NSEAM  has  no  I  and  J  directions  and  so  the  function  value  F  at  1+1  may  not  lie 
next  to  the  function  value  F  at  I.  Thus,  NSEAM  is  written  in  a  very  general  way  which  does  not 
assume  any  direction  and  for  this  reason  it  can  use  any  grid;  this  is  a  significant  advantage  if  one 
expects  to  use  adaptive  or  fully  unstructured  grids. 
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At  this  point,  we  have  begun  modifying  do  loops  inside  of  DIABAT  that  look  like  this: 

Do  J=l, 

Do  1=1,..., IM 

F(I,J)=Some  Operation  (DIRECT ADDRESSING) 

END 

END 

To  loops  that  look  like 
Do  K=l,  ...,NPOIN 

F(K)=Some  Operation  (INDIRECT ADDRESSING) 

END 

Where  IM  and  JM  are  the  number  of  points  in  the  longitude  and  latitude  in  a  lat-lon  grid  and  NPOIN  is 
the  number  of  total  points  in  the  horizontal  (latitude  and  longitude  but  in  an  unstructured  ordering). 

Finally,  to  avoid  too  much  modification  of  both  the  NSEAM  and  DIABAT  codes  it  was  determined 
that  constructing  DIABAT  as  a  separate  library  was  the  best  approach.  This  then  allows  DIABAT  to 
reside  in  different  directories  and  the  linking  of  the  two  codes  is  accomplished  via  library  calls  and 
separate  makefiles. 

IMPACT 

NOGAPS  is  run  operationally  by  FNMOC  and  is  the  heart  of  the  Navy’s  operational  support  to  nearly 
all  DOD  users  worldwide.  This  work  targets  the  next-generation  of  this  system  for  massively  parallel 
computer  architectures.  NSEAM  has  been  designed  specifically  for  these  types  of  computer 
architectures  while  yielding  the  same  high-order  accuracy  as  NOGAPS. 

TRANSITIONS 

Improved  algorithms  for  model  processes  will  be  transitioned  to  6.4  (PE  0603207N)  as  they  are  ready, 
and  will  ultimately  be  transitioned  to  FNMOC  with  future  NOGAPS  upgrades. 

RELATED  PROJECTS 

Some  of  the  technology  developed  for  this  project  will  be  used  immediately  to  improve  the  current 
spectral  transform  formulation  of  NOGAPS  in  other  NRL  proejcts. 
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